﻿using System;
using System.Collections.Generic;
using System.Text;

namespace MLForgeSharp.Models.UnsupervisedLearningModels.DimReduction
{
    /// <summary>
    /// 奇异值分解
    /// </summary>
    public class SVDModel
    {
        private double[,] U; // Left singular vectors
        private double[] S;  // Singular values (diagonal)
        private double[,] Vt; // Right singular vectors (transposed)

        public SVDModel(double[,] matrix, int k)
        {
            // Perform Singular Value Decomposition
            ComputeSVD(matrix, k);
        }

        private void ComputeSVD(double[,] matrix, int k)
        {
            int rows = matrix.GetLength(0);
            int cols = matrix.GetLength(1);

            // Compute A^T * A to get V and S^2
            double[,] ata = MatrixMultiply(Transpose(matrix), matrix);

            // Compute eigenvectors and eigenvalues of A^T * A
            double[] eigenvalues;
            double[,] eigenvectors;
            PowerIteration(ata, k, out eigenvalues, out eigenvectors);

            // Extract singular values and vectors
            S = new double[k];
            Vt = new double[k, cols];
            for (int i = 0; i < k; i++)
            {
                S[i] = Math.Sqrt(eigenvalues[i]);
                for (int j = 0; j < cols; j++)
                {
                    Vt[i, j] = eigenvectors[j, i];
                }
            }

            // Compute U from A * V / S
            U = new double[rows, k];
            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < k; j++)
                {
                    double sum = 0;
                    for (int l = 0; l < cols; l++)
                    {
                        sum += matrix[i, l] * Vt[j, l];
                    }
                    U[i, j] = sum / S[j];
                }
            }
        }

        private void PowerIteration(double[,] matrix, int k, out double[] eigenvalues, out double[,] eigenvectors)
        {
            int n = matrix.GetLength(0);
            eigenvalues = new double[k];
            eigenvectors = new double[n, k];

            Random rand = new Random();
            double[] b = new double[n];
            for (int i = 0; i < n; i++)
            {
                b[i] = rand.NextDouble();
            }

            for (int i = 0; i < k; i++)
            {
                double[] prevB = new double[n];
                double norm = 0;

                // Power iteration to find the dominant eigenvector
                for (int iter = 0; iter < 100; iter++)
                {
                    prevB = (double[])b.Clone();
                    b = MatrixVectorMultiply(matrix, b);

                    // Orthogonalize against previous eigenvectors
                    for (int j = 0; j < i; j++)
                    {
                        double dot = DotProduct(b, GetColumn(eigenvectors, j));
                        for (int l = 0; l < n; l++)
                        {
                            b[l] -= dot * eigenvectors[l, j];
                        }
                    }

                    norm = Norm(b);
                    for (int l = 0; l < n; l++)
                    {
                        b[l] /= norm;
                    }

                    if (Norm(SubtractVectors(b, prevB)) < 1e-10)
                    {
                        break;
                    }
                }

                // Store the eigenvector and eigenvalue
                for (int l = 0; l < n; l++)
                {
                    eigenvectors[l, i] = b[l];
                }
                eigenvalues[i] = norm;

                // Deflate the matrix
                double[,] outerProduct = OuterProduct(b, b);
                for (int l = 0; l < n; l++)
                {
                    for (int m = 0; m < n; m++)
                    {
                        matrix[l, m] -= eigenvalues[i] * outerProduct[l, m];
                    }
                }
            }
        }

        private double[] MatrixVectorMultiply(double[,] matrix, double[] vector)
        {
            int rows = matrix.GetLength(0);
            int cols = matrix.GetLength(1);
            double[] result = new double[rows];

            for (int i = 0; i < rows; i++)
            {
                double sum = 0;
                for (int j = 0; j < cols; j++)
                {
                    sum += matrix[i, j] * vector[j];
                }
                result[i] = sum;
            }

            return result;
        }

        private double DotProduct(double[] a, double[] b)
        {
            double sum = 0;
            for (int i = 0; i < a.Length; i++)
            {
                sum += a[i] * b[i];
            }
            return sum;
        }

        private double[] SubtractVectors(double[] a, double[] b)
        {
            double[] result = new double[a.Length];
            for (int i = 0; i < a.Length; i++)
            {
                result[i] = a[i] - b[i];
            }
            return result;
        }

        private double Norm(double[] vector)
        {
            double sum = 0;
            for (int i = 0; i < vector.Length; i++)
            {
                sum += vector[i] * vector[i];
            }
            return Math.Sqrt(sum);
        }

        private double[] GetColumn(double[,] matrix, int col)
        {
            int rows = matrix.GetLength(0);
            double[] column = new double[rows];
            for (int i = 0; i < rows; i++)
            {
                column[i] = matrix[i, col];
            }
            return column;
        }

        private double[,] OuterProduct(double[] a, double[] b)
        {
            int n = a.Length;
            double[,] result = new double[n, n];
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    result[i, j] = a[i] * b[j];
                }
            }
            return result;
        }

        private double[,] Transpose(double[,] matrix)
        {
            int rows = matrix.GetLength(0);
            int cols = matrix.GetLength(1);
            double[,] result = new double[cols, rows];

            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < cols; j++)
                {
                    result[j, i] = matrix[i, j];
                }
            }

            return result;
        }

        private double[,] MatrixMultiply(double[,] a, double[,] b)
        {
            int rowsA = a.GetLength(0);
            int colsA = a.GetLength(1);
            int colsB = b.GetLength(1);

            double[,] result = new double[rowsA, colsB];

            for (int i = 0; i < rowsA; i++)
            {
                for (int j = 0; j < colsB; j++)
                {
                    double sum = 0;
                    for (int k = 0; k < colsA; k++)
                    {
                        sum += a[i, k] * b[k, j];
                    }
                    result[i, j] = sum;
                }
            }

            return result;
        }

        public double[,] Reconstruct()
        {
            // Reconstruct the original matrix: U * S * Vt
            int rows = U.GetLength(0);
            int cols = Vt.GetLength(1);
            double[,] reconstructed = new double[rows, cols];

            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < cols; j++)
                {
                    double sum = 0;
                    for (int k = 0; k < S.Length; k++)
                    {
                        sum += U[i, k] * S[k] * Vt[k, j];
                    }
                    reconstructed[i, j] = sum;
                }
            }

            return reconstructed;
        }

        public string MatrixToString(double[,] matrix)
        {
            int rows = matrix.GetLength(0);
            int cols = matrix.GetLength(1);
            StringBuilder sb = new StringBuilder();

            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < cols; j++)
                {
                    sb.Append(matrix[i, j].ToString("F2")).Append("\t");
                }
                sb.AppendLine();
            }

            return sb.ToString();
        }
    }

    // 示例程序
    public class SVDModelExample
    {
        public SVDModelExample()
        {
            // Example matrix
            double[,] matrix = 
            {{ 1, 2, 3 },
             { 4, 5, 6 },
             { 7, 8, 9 }};

            // Number of singular values to keep
            int k = 2;

            // Create SVD model
            SVDModel svdModel = new SVDModel(matrix, k);

            // Reconstruct the matrix
            double[,] reconstructedMatrix = svdModel.Reconstruct();

            // Print original matrix
            Console.WriteLine("Original Matrix:");
            Console.WriteLine(svdModel.MatrixToString(matrix));

            // Print reconstructed matrix
            Console.WriteLine("Reconstructed Matrix:");
            Console.WriteLine(svdModel.MatrixToString(reconstructedMatrix));
        }
    }
}
